*********************************************************		
* Replication file for tables and figure in Appendix
* “Reading fiction and economic preferences of rural youth in Burkina Faso”
* Economic Development and Cultural Change 2019
	* Michael Kevane
	* Dept. of Economics
	* Santa Clara University
	* mkevane@scu.edu
*********************************************************		

* Programs might need to install
	*ssc install ingap
	*ssc install estout  

* Programs to find and then install
	*findit matselrc 
	*findit matsave 
	*findit leebounds 
	*findit tsrtest
	*findit permtest2
	*findit leebounds
	*findit strdist

* Set globals 
	clear
	version 12
	set more off 
	dis "`c(username)'"
	* This is for Windows PS; Mac users will have different format
	global user "C:\Users\\`c(username)'\Documents\JTL replication"

* Set language; for tables and figures that alternative titles and labels
	global lang "English"

* Set globals for file structure
	global dta "$user\data"
	global prg "$user\programs"
	global extra "$prg\extra"	
	global pprs "$user\results"
	cd "$user\scrap"
	set more off
	eststo clear

* Create program for randomization inference that calculates means for two groups treatment and control
	capture program drop ttk
	program define ttk, rclass
		args y group
		summ `y' if `group'==1, meanonly
		local fgn=r(mean)
		summ `y' if `group'==0, meanonly
		local dom=r(mean)
		return scalar diff = `fgn'-`dom'
	end

* Loop over two sets of variables: measures of reading, and measures of outcomes of experimental games
	forval z=1(1)2 {
	if "`z'"=="1" {
	* Table of summary statistics BY treatment, with ttest for diff in means, unweighted (not reported in paper)
	loc outcom_jeu livres30 titreyes_mai13 titreyes_aout13 titreyes_mars14 titreyes_mai14 ///
		titrejtlyesnk_mai13 titrejtlyesnk_aout13 titrejtlyesnk_mars14 titrejtlyesnk_mai14   ///
		persjtlyes_aout13 persyes_mars14 persjtlyes_mai14 ///
		faityes_aout13 faityes_mai14
	}
	else if "`z'"=="2" {
	loc outcom_jeu contribconf_mai13 contribconf_aout13 contribconf_mai14  ///
		weighshare_mai13 weighshare_aout13  weighshare_mai14 ///
		contribpub_mai13 contribpub_aout13 contribpub_mai14 ///
		parichoix_mai13 parichoix_aout13 parichoix_mai14 ///
		patiencecount_mai13 patiencecount_aout13 patiencecount_mai14
	}	
	
*********************************************************		
* Table of differences in means - a variety of methods
*********************************************************		
	use  "$dta\JTL_for_EDCC", clear
	eststo clear
	do "$extra\make table of means and ttest results.do" ///  /*this computes means/prop for variables according to dummy 0-1 variables*/
	"treatment"	///  /*these are the categorical variables, the BY variable*/
	" " /// /*these are the 0-1 variables, test for prop diff*/
	"`outcom_jeu' " /// /*these are the continuous variables, ttest for mean diff*/
	"$pprs/App_Tab_temp1" /*file location and name of Excel output*/
	rename v1 varname
	drop if varname==""
	replace varname="method" if varname=="Notes: * p<0.10 ** p<0.05 *** p<0.01."
	rename vd2treatment signif1
	replace signif1="straight averages" if signif1=="" & varname=="method"
	save mean_treat_1, replace

* Table of summary statistics BY treatment, weighted by survey probability, with ttest for diff in means
	use  "$dta\JTL_for_EDCC", clear
	foreach var in `outcom_jeu' {
	replace `var'=`var'*survweight
	}
	eststo clear
	do "$extra\make table of means and ttest results.do" ///  /*this computes means/prop for variables according to dummy 0-1 variables*/
	"treatment"	///  /*these are the categorical variables, the BY variable*/
	" " /// /*these are the 0-1 variables, test for prop diff*/
	"`outcom_jeu' " /// /*these are the continuous variables, ttest for mean diff*/
	"$pprs/App_Tab_temp2" /*file location and name of Excel output*/
	rename v1 varname
	drop if varname==""
	replace varname="method" if varname=="Notes: * p<0.10 ** p<0.05 *** p<0.01."
	rename vd2treatment signif1_w
	rename v2treatment v2treatment_w
	rename v5treatment v5treatment_w
	replace signif1_w="weighted averages" if signif1_w=="" & varname=="method"
	save mean_treat_1w, replace
	
*********************************************************	
* Table of treatment effects where use Stata svy to adjust by survey probabilities
*********************************************************			
	loc catlist treatment
	use  "$dta\JTL_for_EDCC", clear	
	foreach cat in `catlist' {
	foreach var in `outcom_jeu' {
	* calculate means
	quietly svy: quietly mean `var', over(`cat')
	matrix c_`var' = r(table)
	*matrix list c_`var'
	matselrc c_`var' cb_`var' ,  row(b)
	matselrc c_`var' cse_`var',  row(se)
	* do the test that means are equal, test handles svy
	quietly test [`var']0 = [`var']1
	* return list	
	mat U = r(p)*J(1,1,1)
	matrix c1_`var' = (cb_`var',cse_`var',U)
	}
	
	* Now create a list of the variables with a slash so that can append all the matrices
	loc jj=""
	foreach var in `outcom_jeu' {
	loc jj="`jj'  c1_`var'\"
	}
	* have to add one at the end because cannot end on a slash
	loc jj="`jj'  c1_faityes_mai14"
	display "`jj'"
	* use the list to append the matrices
	matrix D = `jj' 
	matrix list D
	
	* Create variables and Stata dataset from matrix
	preserve
	svmat double D	
	forval i=1(1)5 {
	replace D`i'=round(D`i',.001)
	format D`i' %9.3f 
	}
	gen vd=""
	replace vd="*" if D5<=.10
	replace vd="**" if D5<=.05
	replace vd="***" if D5<=.01
	gen E1=""
	loc i=1
	foreach var in `outcom_jeu' {
	replace E1="`var'" if _n==`i'
	loc i=`i'+1
	}
	keep D1 D2 D5 E1 vd
	* D3 and D$ are the standard errors
	drop if D1==.
	drop if E1==""
	gen categ=""
	order categ, before(D1)
	foreach var in D1 D2 D5 E1 vd categ {
	rename `var' `var'_`cat'
	}
	gen E1=E1_`cat'
	save temp_`cat'_means, replace
	restore	
	}
	* merge the various Stata datasets
	loc i=1
	foreach cat in `catlist' {
	if "`i'"=="1" {
	use temp_`cat'_means, clear
	}
	else if "`i"~="1" {
	merge 1:1 E1 using temp_`cat'_means
	drop _merge
	}
	loc i=`i'+1
	}
	order E1, first
	save svymean_est, replace
	
	rename E1 varname
	drop categ_treatment E1_treatment
	drop if varname==""
	rename vd_treatment signif4
	* add an observation at bottom
	  local new = _N + 1
      set obs `new'
	replace varname="method" if varname==""
	replace signif4="survey weighted" if signif4=="" & varname=="method"
	save mean_treat_2, replace
	
	foreach cat in `catlist' {
	erase temp_`cat'_means.dta
	}

*********************************************************	
* Table of treatment effects where use randomization inference for pval 
*********************************************************	
* Use program tsrtest to use randomization inference to see if difference in means significant

	use  "$dta\JTL_for_EDCC", clear	
	gen tsrtestpval=.
	gen ttestpval=.
	gen testvar=""
	gen vd=""
	
	loc i=1
	foreach x in `outcom_jeu' {
	quietly tsrtest treatment r(diff): ttk	`x' treatment
	replace tsrtestpval=r(twotail) if _n==`i' 
	quietly ttest `x', by(treatment)
	replace ttestpval=r(p) if _n==`i'
	replace testvar="`x'" if _n==`i'
	loc i = `i' + 1
	}
	replace vd="*" if tsrtestpval<=.10
	replace vd="**" if tsrtestpval<=.05
	replace vd="***" if tsrtestpval<=.01
	keep tsrtestpval ttestpval testvar vd
	order testvar, first
	list tsrtestpval ttestpval testvar vd in 1/10
	rename testvar varname
	drop if varname==""
	rename vd signif2
	* add an observation at bottom
	  local new = _N + 1
      set obs `new'
	replace varname="method" if varname==""
	replace signif2="RI standard errors" if signif2=="" & varname=="method"
	save mean_treat_3, replace
/*
* Use program tsrtest to use randomization inference to see if difference in means significant, by cluster
* This is muted here because take a LONG time (1-2 days running on ordinary PC desktop in 2018)
	use  "$dta\JTL_for_EDCC", clear
	replace cluster="Bérébaalentours" if cluster=="Béréba alentours"
	gen tsrtestpval=.
	gen ttestpval=.
	gen testvar=""
	gen vd=""
	gen ve=.
	levelsof cluster, local(cl_levels)
	foreach clus in `cl_levels' {
	loc i=1
	foreach x in `outcom_jeu' {
	quietly capture tsrtest treatment r(diff) if cluster=="`clus'": ttk	`x' treatment
	replace tsrtestpval=r(twotail) if _n==`i' 
 	quietly capture ttest `x' if cluster=="`clus'", by(treatment)
	replace ve=r(mu_2)-r(mu_1) if tsrtestpval<=.10 & _n==`i'
	replace ttestpval=r(p) if _n==`i'
	replace testvar="`x'" if _n==`i'
	loc i = `i' + 1
	}
	replace vd="*" if tsrtestpval<=.10
	replace vd="**" if tsrtestpval<=.05
	replace vd="***" if tsrtestpval<=.01
	foreach varx in tsrtestpval ttestpval testvar vd ve {
	gen `varx'_`clus'=`varx' 
	}
	foreach varx in tsrtestpval ttestpval ve {
	replace `varx'=. 
	}
	foreach varx in testvar vd  {
	replace `varx'="" 
	}
	}
	levelsof cluster, local(cl_levels)
	foreach clus in `cl_levels' {
	replace ve_`clus'=round(ve_`clus',.02)
	replace ve_`clus'=round(ve_`clus') if index(testvar_`clus',"contrib")==1
	gen ve_`clus'_str = string(ve_`clus')
	gen ve_`clus'_strsig = ve_`clus'_str+vd_`clus'
	replace ve_`clus'_strsig = "" if ve_`clus'_strsig=="."
	}

	loc cluslist1=""
	foreach clus in `cl_levels' {
	loc cluslist1="`cluslist1' vd_`clus'"
	}
	loc cluslist2=""
	foreach clus in `cl_levels' {
	loc cluslist2="`cluslist2' ve_`clus'"
	}
	loc cluslist3=""
	foreach clus in `cl_levels' {
	loc cluslist3="`cluslist3' ve_`clus'_strsig"
	}	
	keep testvar_Bekuy `cluslist3'
	rename testvar_Bekuy varname
	foreach clus in `cl_levels' {
	rename ve_`clus'_strsig `clus'
	}
	drop if varname==""
	* add an observation at bottom
	  local new = _N + 1
      set obs `new'
	replace varname="method" if varname==""
	save mean_treat_4, replace
*/

**********************************************************************
* Calculate Lee bounds of treatment effect	 
**********************************************************************	
	* Now run Lee bounds with three different tightening variables and all the outcome variables
	* variable bwaba has highest tstat in explaining August 2013 according to regression above
	set more off	 
	loc tightlist  bwaba 
	use  "$dta\JTL_for_EDCC", clear	
	foreach cat in `tightlist' {
	foreach var in  `outcom_jeu' {
	quietly leebounds `var' treatment [pw=survweight], tight(`cat') cieffect vce(bootstrap, reps(250) nodots)
	matrix c_`var' = r(table)
	matselrc c_`var' cb_`var',  row(b) col(treatment:upper) // estimate of coefficient //
	matselrc c_`var' cse_`var',  row(se) col(treatment:upper) // estimate of SE of coefficient //
	matselrc c_`var' cpvalu_`var',  row(pvalue) col(treatment:upper) // estimate of pval of coefficient //
	*matrix list c_`var'
	mat U = e(cilower)*J(1,1,1) // CI for estimated effect //
	mat V = e(ciupper)*J(1,1,1) // CI for estimated effect //
	matrix c1_`var' = (cb_`var',cse_`var',cpvalu_`var',U,V)
	matrix list c1_`var'
	}
	
	* Now create a list of the variables with a slash so that can append all the matrices
	loc jj=""
	foreach var in `outcom_jeu' {
	loc jj="`jj'  c1_`var'\"
	}
	* have to add one at the end because cannot end on a slash
	loc jj="`jj'  c1_faityes_mai14"
	display "`jj'"
	
	* use the list to append the matrices
	matrix D = `jj' 
	matrix list D
	
	* Create variables from matrix
	preserve
	svmat double D	
	forval i=1(1)5 {
	replace D`i'=round(D`i',.001)
	format D`i' %9.3f 
	}
	
	gen vd=""
	replace vd="*" if D3<=.10
	replace vd="**" if D3<=.05
	replace vd="***" if D3<=.01
	
	gen E1=""
	loc i=1
	foreach var in `outcom_jeu' {
	replace E1="`var'" if _n==`i'
	loc i=`i'+1
	}
	
	keep D1 D2 D3 D4 D5 E1 vd
	* D3 and D$ are the standard errors
	drop if D1==.
	drop if E1==""
	gen categ=""
	order categ, before(D1)
	rename D1 coeff
	rename D2 coeff_se
	rename D3 coeff_pval
	rename D4 CI_lower
	rename D5 CI_upperr
	foreach var in coeff coeff_se coeff_pval CI_lower CI_upperr E1 vd categ {
	rename `var' `var'_`cat'
	}
	gen E1=E1_`cat'
	save temp_`cat'_means, replace
	restore	
	}
	
	loc i=1
	foreach cat in `tightlist' {
	if "`i'"=="1" {
	use temp_`cat'_means, clear
	}
	else if "`i"~="1" {
	merge 1:1 E1 using temp_`cat'_means
	drop _merge
	}
	loc i=`i'+1
	}
	order E1, first
	save leebounds_est, replace
	
	foreach cat in `tightlist' {
	erase temp_`cat'_means.dta
	}

	rename E1 varname
	drop if varname==""
	* add an observation at bottom
	  local new = _N + 1
      set obs `new'
	replace varname="method" if varname==""
	*replace signif4="Lee bounds" if signif4=="" & varname=="method"
	save mean_treat_5, replace

*********************************************************	
* Put all the tables together
*********************************************************	
	use mean_treat_5, clear
	merge 1:1 varname using  mean_treat_1w
	drop _merge
	merge 1:1 varname using   mean_treat_2
	drop _merge
	merge 1:1 varname using   mean_treat_3
	drop _merge	
	*merge 1:1 varname using   mean_treat_4
	*drop _merge	
	merge 1:1 varname using   mean_treat_5
	drop _merge
	gen ordertab=.
	if "`z'"=="1" {
	replace ordertab=1 if varname=="livres30"
	replace ordertab=2 if varname=="atleast1liv30"
	replace ordertab=3 if varname=="atleast1titrys"
	replace ordertab=4 if varname=="titreyes_mai13"
	replace ordertab=5 if varname=="titreyes_aout13"
	replace ordertab=6 if varname=="titreyes_mars14"
	replace ordertab=7 if varname=="titreyes_mai14"
	replace ordertab=8 if varname=="titrejtlyesnk_mai13"
	replace ordertab=9 if varname=="titrejtlyesnk_aout13"
	replace ordertab=10 if varname=="titrejtlyesnk_mars14"
	replace ordertab=11 if varname=="titrejtlyesnk_mai14"
	replace ordertab=12 if varname=="persyes_aout13"
	replace ordertab=13 if varname=="persyes_mars14"
	replace ordertab=14 if varname=="persyes_mai14"
	replace ordertab=15 if varname=="faityes_aout13"
	replace ordertab=16 if varname=="faityes_mars14"
	replace ordertab=17 if varname=="faityes_mai14"
	replace ordertab=18 if varname=="readingpercorr_aout13"
	replace ordertab=19 if varname=="readingpercorr_mai14"
	replace ordertab=20 if varname=="treatment"
	replace ordertab=21 if varname=="Sample size"
	replace ordertab=23 if varname=="method"
	}
	else if "`z'"=="2" {
	replace ordertab=1 if varname=="contribconf_mai13"
	replace ordertab=2 if varname=="contribconf_aout13"
	replace ordertab=3 if varname=="contribconf_mai14"
	replace ordertab=4 if varname=="weighshare_mai13"
	replace ordertab=5 if varname=="weighshare_aout13"
	replace ordertab=6 if varname=="weighshare_mai14"
	replace ordertab=7 if varname=="contribpub_mai13"
	replace ordertab=8 if varname=="contribpub_aout13"
	replace ordertab=9 if varname=="contribpub_mai14"
	replace ordertab=10 if varname=="parichoix_mai13"
	replace ordertab=11 if varname=="parichoix_aout13"
	replace ordertab=12 if varname=="parichoix_mai14"
	replace ordertab=13 if varname=="patiencecount_mai13"
	replace ordertab=14 if varname=="patiencecount_aout13"
	replace ordertab=15 if varname=="patiencecount_mai14"
	}
	sort ordertab
	export excel using "$pprs/TabRILee`z'", sheetreplace firstrow(var) nolabel

	erase "mean_treat_1w.dta"
	erase "leebounds_est.dta"
	erase "svymean_est.dta"
	forval i=1(1)3 {
	erase "mean_treat_`i'.dta"
	}
	}
